Circadian-period variation underlies the local adaptation of photoperiodism in the short-day plant Lemna aequinoctialis

Summary Phenotypic variation is the basis for trait adaptation via evolutionary selection. However, the driving forces behind quantitative trait variations remain unclear owing to their complexity at the molecular level. This study focused on the natural variation of the free-running period (FRP) of the circadian clock because FRP is a determining factor of the phase phenotype of clock-dependent physiology. Lemna aequinoctialis in Japan is a paddy field duckweed that exhibits a latitudinal cline of critical day length (CDL) for short-day flowering. We collected 72 strains of L. aequinoctialis and found a significant correlation between FRPs and locally adaptive CDLs, confirming that variation in the FRP-dependent phase phenotype underlies photoperiodic adaptation. Diel transcriptome analysis revealed that the induction timing of an FT gene is key to connecting the clock phase to photoperiodism at the molecular level. This study highlights the importance of FRP as a variation resource for evolutionary adaptation.


INTRODUCTION
Understanding the molecular framework underlying the phenotypic changes in local adaptation is an important goal in evolutionary ecology. Natural variation in phenotypic traits plays a key role in evolutionary adaptation (Barret and Schluter, 2008;Alonso-Blanco et al., 2009;Matuszewski et al., 2015). However, the molecular mechanisms that generate such phenotypic variations in natural settings are still poorly understood because they are complex traits at the molecular level (Boyle et al., 2017;Sella and Barton, 2019;Feiner et al., 2021).
The circadian clock is an endogenous timing system that generates a self-sustained daily oscillation entrained to day-night cycles via environmental cues and regulates various physiologies to optimize the timing of development, stress responses, and metabolism (Pittendrigh 1960;Greenham and McClung, 2015). The eukaryotic circadian clock comprises feedback loops of multiple clock genes (Saini et al., 2019). The free-running period (FRP) of the circadian clock is a polygenic trait that exhibits natural variations (Michael et al., 2007;Brown et al., 2008;Pivarciova et al., 2016;Salmela and Weinig, 2019). The clock with an FRP deviating from 24 h is entrained to the 24-h day-night cycles (Johnson et al., 2003). The peak timing (phase) of an entrained rhythm also shows a natural variation. Previous studies suggest that diel phase determination is related to FRP in various organisms; a longer FRP results in a phase delay phenotype (Aschoff and Pohl, 1978;Ré mi et al., 2010;Graf et al., 2010;Granada et al., 2013). Because the clock phase is relayed to physiological outputs, we hypothesized that FRP variation is associated with the local adaptation of the phase phenotype of clock-dependent physiology (Lankinen 1993;Dominoni et al., 2013;Helm et al., 2017).
FRP variation within 20-28 h has been reported in several plant species. In Arabidopsis thaliana, FRP showed a latitudinal cline and quantitative trait loci, including several PSEUDO-RESPONSE REGULATOR (PRR) genes, which encode a core component of the circadian clock, were detected (Michael et al., 2003). In Swedish A. thaliana accessions, an allele of the COLD REGULATED 28 (COR28) gene of accessions in southern areas was associated with long-period and late-flowering phenotypes (Rees et al., 2021). Latitudinal clines of FRP in wild monkey flower (Mimulus guttatus) and cultivated soybean (Glycine max) were also reported (Greenham et al., 2017). During the domestication of tomato, an allele of the EMPFINDLICHER IM DUNKELROTEN LICHT 1 (EID1) gene for a longer FRP was selected because of its higher crop performance under long-day conditions (Mü ller et al., 2016). These results imply that FRP variation is linked to local adaptation of clock-dependent temporal traits. However, direct evidence of the association between FRP and adaptive temporal traits is still limited because of difficulties in evaluating the latter in natural populations.
As a clock-dependent physiology in plants, photoperiodic flowering is a typical target of local adaptation; the critical day length (CDL) often shows a latitudinal cline to adjust the flowering time to local habitat conditions (Ray and Alexander 1966;Imamura et al., 1966;Katayama 1977;Yukawa and Takimoto, 1976). In Arabidopsis, the CONSTANS (CO) gene plays a key role in photoperiodism (Suá rez-Ló pez et al., 2001). According to the external coincidence model, its diel expression pattern is critical in CDL determination, and the coincidence of light with CO expression activates the FLOWERING LOCUS T (FT) gene, the master regulator of floral induction. A clock-related gene, GIGANTEA (GI), is responsible for the expression pattern of CO (Sawa et al., 2007). Notably, an Arabidopsis clock gene mutant with a shorter FRP exhibited an early flowering phenotype under short-day conditions. A shorter FRP resulted in earlier CO induction, which broadened the coincidence window to induce FT expression earlier (Yanovsky and Kay, 2002). Thus, the induction timing of FT is a potential target for the local adaptation of CDL. Consistent with this idea, Populus trichocarpa showed a latitudinal cline in flowering and growth cessation, and the induction timing of PtFT1 was earlier in the southern populations (Bö hlenius et al., 2006). Because the natural variation of FRP is expected to be related to clock phase variation, we hypothesize that FPR-dependent phase changes in photoperiodic genes alter the coincidence window for FT activation and contribute to the local adaptation of CDL.
To investigate the relationship between FRP and CDL, we focused on the short-day duckweed Lemna aequinoctialis. Duckweed species are small, free-floating aquatic plants that have been studied in diverse fields ranging from ecology to genomics (Acosta et al., 2021). High-throughput phenotyping using the particle bombardment transfection technique and a luciferase reporter can be used in duckweed species to study the variation in circadian traits (Miwa et al., 2006;Muranaka et al., 2015;Isoda et al., 2022). A latitudinal cline of the CDL for photoperiodic flowering has been described in L. aequinoctialis Japanese populations (Yukawa and Takimoto, 1976). L. aequinoctialis produces flower buds within one week of photoperiodic treatment, allowing rapid CDL determination. Isolated duckweed strains are maintained by clonal growth, and laboratory experiments using clonal strains have enabled the precise determination of multiple phenotypic traits of the original genotypes from natural populations. By taking advantage of this plant material, we detected a significant correlation between FRPs and locally adaptive CDLs.

Strain collection from natural populations in Japan
We collected 72 strains of L. aequinoctialis from 20 populations between latitudes 31.5 N and 43.8 N (Figure 1A; Table S1). To estimate both intra-and inter-population variations, 11-12 strains were collected from four populations (N32Ka, N35Ht, N35So, and N44Ha). The strains were aseptically maintained by clonal growth. The frond size showed intra-and inter-population variation, suggesting heterogeneous genotypes, even within a population ( Figure S1). The estimated genome sizes of all strains were similar ( Figure S2).

Natural variation in CDLs
These strains detected day lengths with a resolution of 0.5 h, and their CDLs ranged from 11.2 to 14.2 h ( Figures 1B, 1C, and 1D). As previously reported (Yukawa and Takimoto, 1976), a latitudinal cline was observed in CDLs: 11.2-12.6 h in the south, 12.4-14.2 h in the north, and a large variation in CDL was observed at approximately 35 N ( Figure 1E). The longer CDL at higher latitudes is consistent with that of other short-day flowering species (Ray and Alexander 1966;Imamura et al., 1966;Katayama 1977). Intrapopulation differences of 1-1.5 h in CDLs were observed in each population of N32Ka, N35Ht, N35So, and N44Ha ( Figure 1F). The significant differences between these populations suggest a geographical differentiation of CDLs. The difference in CDL between the two 35 N populations, N35Ht and N35So, suggests that the adaptive CDL at each site may differ, even at the same latitude ( Figure 1F). The flooding period of the paddy fields differed between the sampling sites depending on the cultivation schedule of the different rice cultivars. The longer CDL of N35So, corresponding to earlier flowering, is consistent with the earlier cessation of flooding at this site than at the N35Ht site ( Figure 1G). The CDL of the N32Ka population was shorter than that of the N35 population ( Figure 1F). The CDLs of the three populations appear to be linked to the cultivation schedule at each site ( Figure 1G). The complete drainage of paddy fields during ll OPEN ACCESS  Figure 2A). All strains showed clear diel luminescence rhythms with morning peaks at 15L9D. Among the four populations, the peak time of N32Ka was significantly later than those of the other three populations ( Figure 2B). To estimate FRP under LL, a fast Fourier transform-nonlinear least squares (FFT-NLLS) analysis was used (Zielinski et al., 2014), in which the rhythm significance was estimated by a relative amplitude error (RAE) that increased from 0 to 1 as the rhythm neared statistical insignificance ( Figure 2C). Seven strains with a high RAE value (>0.1) or a high SD of FRP (>1.5 h) were excluded from the subsequent analysis ( Figure S3). The FRP of N32Ka was significantly longer than those of the other populations, corresponding to its late peak time ( Figures 2B and 2D). This is consistent with the idea of an FRPdependent phase phenotype, that is, a positive correlation between FRP and peak time (Helm et al., 2017).

A negative correlation between FRP and CDL
Correlation analysis using data from all 72 strains (20 populations) suggested an effect of FRP-dependent phase phenotype on photoperiodism ( Figure 3A). A positive correlation between FRP and peak time was observed at the genotypic level. In addition, CDL showed a negative correlation with both FRP and peak time. These results suggest that the variation in the FRP-dependent phase phenotype underlies the  Figure 3B). In this model, the FT induction is assumed to be gated by the circadian clock and permitted in the dark. Thus, FT expression responds to night lengthening. A longer FRP causes a phase delay in the gate timing and, consequently, a shorter CDL, resulting in a negative correlation between FRP and CDL.
Diel expression patterns of FT, CO, and GI homologs We obtained diel transcriptomes of L. aequinoctialis pure line (Nd strain, Muranaka et al., 2015) under three photoperiod conditions to explore FT homologs. The de novo assembly detected five FT homologs. They were named LaFTh1-5 ( Figures S4A and S4B). In a phylogenetic tree with FT homologs of various plant species, LaFTh1 composes a clade with a Zea mays FT homolog (ZmZCN8) and a Brachypodium distachyon FT homolog (BdFTL9). LaFTh1 was induced before the night's end, and its expression increased during longer-night conditions ( Figures 3C and S4C). Interestingly, it has been reported that ZmZCN8 and BdFTL9 are induced during the night only under short-day conditions (Meng et al., 2011;Qin et al., 2019). The transgenic A. thaliana carrying an LaFTh1 overexpression construct showed an early flowering phenotype, suggesting its floral induction activity ( Figure S4D). LaFTh1 was induced even in the extended dark following a non-flowering long-day condition, suggesting that LaFTh1 expression directly responds to night lengthening. The timing of LaFTh1 induction differed among the three strains with different CDLs, and the variation in induction timing appeared to be responsible for the CDL variation ( Figures 3D and 3E). Thus, the iScience Article induction timing of LaFTh1 is key for connecting the clock phase to photoperiodism at the molecular level. Notably, the LaFTh1 expression levels showed a large difference. LaFTh1 expression in the N44Ha08 strain was 10-fold higher than in the other two strains ( Figure 3E). This high basal level may have caused the extremely long CDL of 14.2 h ( Figure 1F).
Next, we explored the L. aequinoctialis homologs of CONSTANS (CO) and GIGANTEA (GI ), which are core components of day-length measurement, by regulating FT expression in Arabidopsis photoperiodism (Sawa et al., 2007). We detected three CO homologs (LaCOh1-3) with amino acid sequence similarities to AtCOL3 and AtCOL4, rather than AtCO ( Figure S5A). The LaCOh1-3 genes showed diel expression patterns, but none had the midnight peak phase observed for AtCO and OsHd1 (the rice CO ortholog) (Yanovsky and Kay 2002;Hayama et al., 2003). In contrast, a GI homolog (LaGIh1) showed diel expression patterns that were similar to those of Arabidopsis and rice ( Figure S5B) (Sawa et al., 2007;Hayama et al., 2003). These results suggest that the molecular mechanism of day-length measurement in L. aequinoctialis is different from that in Arabidopsis and rice.

DISCUSSION
Our study demonstrated the effect of the FRP-dependent phase phenotype on locally adapted traits; CDL variation is tightly coupled with FRP variation, even in natural populations. Diel transcriptome analysis suggests that FT induction timing, which is likely to be related to FRP-dependent phase variation, is a determining factor for CDL and, consequently, the flowering season. These results indicated the importance of circadian phenotypes in the natural selection of seasonal phenotypes. This study also provides new insights into the local adaptation of CDLs to paddy field environments. Consistent with previous studies using short-day plants, the populations of L. aequinoctialis in Japan showed a latitudinal cline in their CDLs, suggesting that the flowering season of this plant adapted to the local climate, exhibiting a latitudinal cline in temperature (Hut et al., 2013). In addition to the latitudinal cline, another factor in the local adaptation of CDLs should be considered to explain the variation of CDLs among populations at similar latitudes (Figure 1E). Interestingly, the CDLs appeared to adapt to the flooding season at each sampling site ( Figure 1G). The improved drainage systems associated with the agricultural modernization of Japan that began in the 1950s have resulted in dry paddy field environments in winter and consequently altered the weed flora (Katayama et al., 2015). The natural variation of the CDLs possibly contributed to the adaptive evolution of the flowering season in response to the rapid artificial shift of the flooding season. In addition to interpopulation variations, CDLs showed intrapopulation variations ( Figure 1F). The flooding period of paddy fields is highly dependent on the rice cultivar used. Rice cultivars planted in the same area may be spatially or yearly heterogeneous. Such artificial fluctuations may promote intrapopulation variation in CDLs. The natural variation in FRP-dependent phase phenotypes may play a role in maintaining CDL variations. In A. thaliana, the recombinant inbred lines of two accessions with similar FPRs showed variation as substantial as that observed in the global collection (22.0-28.5 h) (Michael et al., 2003). Such transgressive segregation of FPR is caused by the complex interference among many quantitative loci (Rieseberg et al., 2003), and may play an important role in providing phenotypic variations for rapid evolutionary adaptation to irregular environmental changes (Hamann et al., 2021). In this respect, the apparent tetraploidy of L aequinoctialis ( Figure S2; Beppu and Takimoto, 1981;Wang et al., 2011) could have contributed to the generation of phenotypic variation by increasing the rates of evolution and number of components in gene regulatory networks (Sé mon and Wolfe 2007; Selmecki et al., 2015).
As the circadian clock regulates various aspects of physiology (Greenham and McClung, 2015;Helm et al., 2017), FRP variation potentially contributes to the fine-tuning of many different traits in the process of local adaptation. Conversely, selection pressure on clock-dependent physiological traits may favor mutations in clock-related genes, which alter FRP, and consequently, the diel phase (Figure 4). This is a possible reason for the maintenance of FRP variation with standing genetic variations in clock-related genes in the natural populations (Hut et al., 2013;Salmela and Weinig, 2019). Thus, the adaptive significance of FRP variation should be considered in relation to the phase phenotype under day-night conditions. Interestingly, a latitudinal cline of FRP has been reported in the linden bug Pyrrhocoris apterus, in which circadian clock genes are responsible for the photoperiodic response of the reproductive diapause (Pivarciova et al., 2016;Urbanová et al., 2016). FRP-dependent phase variation may be involved in local adaptation in various organisms. Longer FRP in soybean and tomato appears to be associated with higher crop performance under long-day conditions (Mü ller et al., 2016;Greenham et al., 2017 iScience Article day-night cycles and has evolved as a hub of environmental response systems involving multiple input pathways (Dixon et al., 2014). As a result, the circadian system can function as a source of variation for temporal traits, resulting from many quantitative loci for FRP in the process of adaptive evolution (Nagel and Kay, 2012;Boyle et al., 2017;Sun et al., 2021).

Limitations of the study
Although L. aequinoctialis is a good model for investigating the association between FRP and CDL, the understanding of the molecular basis of photoperiodism in this plant is insufficient to directly verify the effect of the FRP-dependent phase phenotype on CDL determination at the molecular level. To access the genetic basis of the natural variation of FRP and CDL, functional analysis of homologs of the circadian clock, light signaling, and photoperiodism genes is required. We should note the importance of factors other than FRP in CDL determination because the high basal level of LaFTh1 appears to contribute to the longer CDL in the northern strain. Future studies should provide genomic information and genetic tools for L. aequinoctialis. In addition to molecular analysis, an ecological field study is required to confirm CDL adaptation to the flooding period in paddy fields.

Materials availability
This study did not generate new unique materials or reagents.

Plant materials and growth conditions
Seventy-two L. aequinoctialis strains were collected from 20 sites in Japan (Table S1). The collected plants were sterilized by washing with 5% sodium hypochlorite for several minutes, followed by washing with sterilized water. From the successfully sterilized plants, a single colony was isolated. Only one strain was isolated from a paddy field, and a maximum of 12 paddy fields were selected at each site. These strains were maintained aseptically in a growth medium (NF medium containing 1% sucrose) under long-day conditions (15L9D; approximately 50 mmol m À2 s À1 ) in an incubator (LH350-SP, NK system, Japan) at 25 G 1 C  L. aequinoctialis 6746 and Lemna gibba p8L strains were maintained in the same conditions. Before the experiments, the plants were grown in a flask or 6-well plastic plate filled with 30 mL or 8 mL of growth medium, respectively.

Genome size estimation
The relative genome sizes of L. aequinoctialis strains were estimated based on 4',6-diamidino-2-phenylindole diacetate dye flow cytometry (DAPI-FCM) using a Partec CyStain UV Precise P reagent kit (Sysmex iScience Article Partec GmbH, Germany) and a flow cytometer (Partec Ploidy Analyzer PA-I, Sysmex Partec GmbH). Lemna gibba p8L was used as the standard. The relative genome sizes of the L. aequinoctialis strains were calculated as the relative peak positions. The genome size of L. gibba has been estimated to be 486 Mbp (Wang et al., 2011).

Frond-length measurements
The top view image of the plants was captured using a digital camera (EOS 5D mark3 [Canon, Japan] with SP 90 mm Di MACRO [Tamron, Japan] or EOS Kiss X5 with EF-S55-250 mm [Canon]). The frond length was manually measured as the length from the base to the apex of the mother frond using ImageJ 1.53c software ( Figure S1A). The frond lengths of the eight colonies were measured for each strain.

Measurement of the critical day length in photoperiodic flowering
A colony consisting of three or four visible fronds was placed in each well of a 24-well plastic plate filled with 2 mL of growth medium. These plants were subjected to a one-week photoperiodic treatment with day lengths of 9, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5, and 15 h in an eight-chamber incubator (LH30-8-CT, NK system, Japan). The temperature was maintained at 25 G 1 C and the light intensity was approximately 30 mmol m À2 s À1 . After treatment, the plants were bleached with 70% ethanol. The number of flowering and non-flowering fronds was counted for each colony under a stereoscopic microscope. Flowering rate was calculated as the percentage of flowering fronds to the total number of fronds in at least two replicates. The critical day length was determined as the day length where 50% of the maximum flowering rate (Fmax) was expected in a piecewise linear function for the obtained flowering rate ( Figure 1D).

Meteorological data
Temperature data at each site were obtained from the Automated Meteorological Data Acquisition System (AMeDAS) data download service (https://www.data.jma.go.jp/risk/obsdl/index.php). Day lengths were calculated using the online service of the National Astronomical Observatory of Japan (https://eco.mtk. nao.ac.jp/cgi-bin/koyomi/koyomix_en.cgi).

Monitoring luminescence rhythms
Luminescence monitoring using a circadian luminescent reporter was performed as follows . Plasmid DNA carrying the luciferase reporter gene pUC-AtCCA1::LUC+ (AtCCA1::LUC, Nakamichi et al., 2004) was coated on 0.48 mg of gold particles (1.0 mm diameter, Bio-Rad) and introduced into plants laid on a 60-mm plastic dish using a particle bombardment system (PDS-1000/He, Bio-Rad) according to the manufacturer's instructions (vacuum, 27 mmHg; helium pressure, 450 psi). After particle bombardment, the plants were divided into three 35-mm dishes (approximately four colonies per dish) filled with 4 mL of growth medium containing D-luciferin (0.2 mM potassium salt, Wako). A luminescence dish-monitoring system with photomultiplier tubes (H7360-01; Hamamatsu Photonics K.K., Japan) was used for the luminescence measurements. To reduce the background chlorophyll fluorescence, a shortpass filter (SVO630; Asahi Spectra, Japan) was set at the detection site of the photomultiplier tubes. Each dish was subjected to 30 s of measurements every 20 min. The monitoring system was placed in an incubator (KCLP-1000I-CT; NK System, Japan) with fluorescent lamps (FL20SSW/18; Mitsubishi Electric Co., Japan). The temperature was maintained at 25 G 1 C and the light intensity was approximately 30 mmol m À2 s À1 .

Time-series analysis
Peak detection and FRP estimation were performed as follows (Muranaka and Oyama, 2016). The peak positions were estimated by local quadratic curve fitting. For period estimation, the obtained luminescence time series were detrended by subtracting the 24-h moving average. Thereafter, the amplitude was normalized by dividing by the 24-h moving standard deviation. The normalized time series of 60 to 132 h was analyzed using fast Fourier transform-nonlinear least squares (FFT-NLLS) to determine the FRP (Zielinski et al., 2014). FFT-NLLS was based on a multicomponent cosine fit. Rhythm significance was estimated by a relative amplitude error (RAE) that increases from 0 to 1 as the rhythm approached statistical insignificance. These analyses were performed using R scripts developed and run with R 3.6.3 (http:// r-project.org/).